## Non-racial attitudes (whites only)
## Brian T. Hamel and Bryan Wilcox-Archuleta
## First: 24 September 2019 
## Last: 19 March 2020

## Loading packages
## install.packages(c("estimatr", "tidyverse", "magrittr", "texreg", "gridExtra", "scales"))
library(estimatr)
library(tidyverse)
library(magrittr)
library(texreg)
library(gridExtra)
library(scales)

## Loading data
load("01_data/dta.RData")

## Create shell
shell = dta %>% filter(white == 1) %$%
  expand.grid(racial_flux = seq(min(racial_flux, na.rm = TRUE), max(racial_flux, na.rm = TRUE), by = 1),
              pid7 = round(mean(pid7, na.rm = TRUE), digits = 2),
              ideo5 = round(mean(ideo5, na.rm = TRUE), digits = 2),
              female = round(mean(female, na.rm = TRUE), digits = 2),
              age = round(mean(age, na.rm = TRUE), digits = 2),
              faminc = round(mean(faminc, na.rm = TRUE), digit = 2), 
              educ = round(mean(educ, na.rm = TRUE), digits = 2), 
              pct_white = round(mean(pct_white, na.rm = TRUE), digits = 2), 
              pct_black = round(mean(pct_black, na.rm = TRUE), digits = 2),
              pct_unemployed = round(mean(pct_unemployed, na.rm = TRUE), digits = 2),
              pct_college = round(mean(pct_college, na.rm = TRUE), digits = 2),
              log_per_cap_inc = round(mean(log_per_cap_inc, na.rm = TRUE), digits = 2),
              gini = round(mean(gini, na.rm = TRUE), digits = 2),
              south = round(mean(south, na.rm = TRUE), digits = 2),
              non_rural = round(mean(non_rural, na.rm = TRUE), digits = 2),
              log_pop_density = round(mean(log_pop_density, na.rm = TRUE), digits = 2)) %>% 
  na.omit()

## Models and predicted probabilities
abortion = lm_robust(abortion ~ racial_flux + pid7 + ideo5 + female + age + faminc
                     + educ + pct_white + pct_black + pct_unemployed
                     + pct_college + log_per_cap_inc + gini + south + non_rural
                     + log_pop_density, data = dta %>% filter(white == 1),
                     clusters = zipcode, se_type = "stata")

pred_abortion = cbind(predict(abortion, shell, 
                              se.fit = TRUE, type = "response"), 
                      shell)

climate = lm_robust(climate ~ racial_flux + pid7 + ideo5 + female + age + faminc
                    + educ + pct_white + pct_black + pct_unemployed
                    + pct_college + log_per_cap_inc + gini + south + non_rural
                    + log_pop_density, data = dta %>% filter(white == 1),
                    clusters = zipcode, se_type = "stata")

pred_climate = cbind(predict(climate, shell, 
                             se.fit = TRUE, type = "response"), 
                       shell)

guns = lm_robust(guns ~ racial_flux + pid7 + ideo5 + female + age + faminc
                 + educ + pct_white + pct_black + pct_unemployed
                 + pct_college + log_per_cap_inc + gini + south + non_rural
                 + log_pop_density, data = dta %>% filter(white == 1),
                 clusters = zipcode, se_type = "stata")

pred_guns = cbind(predict(guns, shell, 
                          se.fit = TRUE, type = "response"), 
                  shell)

## Table of coefs., and save
##############
## TABLE A5 ##
##############
texreg(list(abortion, climate, guns),
       file = "03_output/non_racial.tex", 
       label = "non_racial",
       caption = "Racial Flux and Non-Racial Attitudes (Whites)", 
       custom.model.names = c("\\textit{Abortion}", "\\textit{Climate Change}", 
                              "\\textit{Gun Control}"),
       custom.coef.names = c("Intercept", "Racial Flux", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE)

## Plot, and save
pred_abortion = cbind(pred_abortion, outcome = "Abortion")
pred_climate = cbind(pred_climate, outcome = "Climate Change")
pred_guns = cbind(pred_guns, outcome = "Gun Control")
pred_att = bind_rows(pred_abortion, pred_climate, pred_guns) %>% 
  mutate(upper = fit + 1.96 * se.fit,
         lower = fit - 1.96 * se.fit)

###############
## FIGURE A1 ##
###############
att_plot = ggplot(pred_att, aes(x = racial_flux, y = fit, ymin = lower, ymax = upper)) +
  geom_line(color = "red4") +
  geom_ribbon(alpha = .2, fill = "red1") +
  facet_wrap(~ outcome, nrow = 2, ncol = 2, scales = "free") +
  labs(y = "Predicted Attitude", 
       x = "Racial Flux") +
  geom_rug(data = dta %>% filter(white == 1), aes(x = racial_flux), inherit.aes = FALSE, sides = "b") +
  scale_y_continuous(labels = number_format(accuracy = 0.01)) +
  theme(legend.title = element_blank(),
        panel.spacing = unit(1, "lines"),
        axis.line.y = element_blank()) +
  ggsave(file = "03_output/non_racial.png", height = 4, width = 4, units = "in", dpi = 600)

## Clear R 
rm(list = ls())